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1.  INTRODUCTION 


The  subject  of  the  flow  behind  a  projectile  in  flight  has  been  studied  extensively  for 
many  years.  Since  the  drag  on  the  projectile  due  to  the  reduced  pressme  on  the  base  is  a 
significant  portion  of  the  total  drag,  aerodynamicists  have  devised  various  methods  for 
reducing  the  "base  drag".  An  important  technique  for  reducing  the  base  drag  (i.e., 
increasing  the  base  pressure)  is  the  injection  of  combustible  gases  from  the  base.  These 
gases  subsequently  mix  with  the  free  stream  air  and  bum  downstream  of  the  projectile. 
This  method  for  reducing  drag  was  first  suggested  many  decades  ago  (e.g.,  Baker,  Davis 
and  Matthews  1951).  A  collection  of  papers  on  analytic  and  experimental  studies  of  base 
combustion  was  edited  by  Murthy  et  al.  (1976).  This  work  also  includes  a  review  of  base 
flow  phenomena  with  and  without  injection  by  Murthy  and  Osborn  (1976)  through  1974. 
Numerous  approximate  techniques  for  analysis  of  the  base  combustion  flow  problem, 
and  the  influence  on  base  drag  were  presented.  Strahle  and  his  co-workers,  (Hubbartt, 
Strahle  and  Neale  1981  and  Strahle,  Hubbartt  and  Walterick  1982),  have  experimentally 
studied  base  burning  and  external  burning  in  supersonic  flow  using  and  diluents. 

The  effect  of  injectant  molecular  weight  and  energy  content  on  base  drag  was 
investigated. 

The  increase  in  capability  for  analyzing  complicated  flow  problems  using 
computational  fluid  dynamics  (CFD)  techniques,  and  the  availability  of  super  computers 
have  led  to  improved  numerical  analysis  of  both  forebody  and  base  flow  problems. 
Sturek,  Nietubicz,  Sahu,  Danberg  and  others  (Sturek  et  al.  1978;  Nietubicz,  Inger  and 
Danberg  1984;  Sahu,  Nietubicz  and  Steger  1985;  Sahu  1986;  and  Sahu  and  Danberg 
1986)  from  the  U.S.  Army  Ballistic  Research  Laboratory  have  utilized 
inviscid/boundaiy-layer  coupled  techniques  and  implicit  Navier-Stokes  codes  (Nietubicz, 
Pulliam  and  Steger  1980)  to  study  the  flow  fields  for  many  different  projectile 
configurations.  These  works  have  considered  base  flows  without  injection  as  well  as  with 
injection  of  cold  or  hot  air.  Sahu  and  Nietubicz  (1984)  and  Childs  and  Caruso  (1987) 
have  also  considered  the  base  flow  problem  with  a  propulsive  jet  However,  the  present 
work  concentrates  on  the  so-called  base  bleed  phenomena  m  which  only  a  relatively 
small  mass  of  gas  is  injected  from  the  base. 

Modem  U.S.  Army  projectiles  utilize  injection  gases  generated  by  burning  a  fuel 
rich  solid  propellant  whose  primary  combustion  products  are  H, ,  CO,  HQ  and  other 
noncombustible  gases.  These  injection  gases  exit  the  projectile  base  at  low  speed 
relative  to  the  initial  flight  speed,  and  the  duration  of  injection  is  of  order  30  seconds. 

No  detailed  analysis  technique  has  been  developed  yet  for  the  base  flow  combustion 
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problem.  The  present  effort  develops  several  combustion  models  suitable  for  inclusion 
in  Navier-Stokes  computational  procedures  for  projectile  base  flow  field  prediction. 
These  models  have  evolved  from  the  hydrogen-air  combustion  literature  for  scramjet 
and  ramjet  reacting  flow  problems,  and  from  the  hydrocarbon  combustion  literature. 

Hydrogen  combustion  has  been  studied  extensively  for  many  years.  For  example, 
Spiegler,  Wolfshtein  and  Manheimer-Timnat  (1976)  utilized  a  seven  species,  eight 
reaction  model  including  the  influence  of  turbulent  fluctuations.  Janicka  and  KoUmann 
(1979)  proposed  a  two-scalar  formulation  based  on  a  seven  reaction  system  and  a 
two-dimensional  pdf  for  modeling  the  effect  of  turbulence  in  an  H,  -air  difhision  flame. 
This  model  assumes  that  the  two-body  "shuffle"  reactions  occur  very  rapidly  so  that  they 
are  in  equilibrium,  while  the  slower  three-bo^  recombination  reactions  are  considered 
kinetically. 

Rogers  and  Chinitz  (1983)  developed  a  two-step  global  reaction  model  for  H^-air 
combustion  at  one  atmosphere  pressure.  This  model  requires  only  five  species  including 
N, ;  therefore,  it  is  more  efficient  than  the  more  extensive  mechanisms.  Also,  this  model 
includes  the  effect  of  stoichiometiy  on  the  global  reaction  rates.  Uenishi,  Rogers  and 
Northam  (1987)  used  the  Rogers  and  Chinitz  model  successfully  for  three-dimensional 
predictions  behind  a  back-step  in  a  supersonic  combustor.  More  recently,  Jachimowski 
(1988)  developed  a  13  species,  33  reaction  model  for  H,  -air  combustion  studies  in 
hypersonic  flows  over  a  range  of  initial  temperatures.  A  nine  species,  18  reaction  model 
was  also  proposed  by  Jachimowski  (1988). 

Evans  and  Schexnayder  (1980)  used  the  Spiegler,  Wolfshtein  and  Manheimer- 
Timnat  (1976)  reaction  system  and  a  12  species,  25  reaction  system  alo  ig  with  the 
"unmixedness"  formulation  of  Spiegler  to  conq)are  with  several  different  supersonic 
flame  test  cases.  The  important  conclusions  from  this  study  were  that  the  25  reaction 
system  was  superior  to  the  eight  reaction  system  for  the  prediction  of  ignition,  but  that 
otherwise  the  eight  reaction  system  was  acceptable.  Unmixedness  also  had  a  significant 
influence  in  one  case  where  ignition  failed  to  occur;  otherwise,  the  effect  was  moderate. 
Eklund,  Drummond  and  Hassan  (1990)  used  a  modified  seven  reaction  set  patterned 
after  that  of  Jachimowski  (1988)  to  calculate  the  combustion  in  turbulent  shear  layers 
and  compare  with  experimental  data. 

The  consideration  of  H,  and  CO  in  flames  has  not  been  as  extensive  as  hydrogen 
alone.  Early  work  was  performed  at  Princeton  University  by  Diyer  (1972)  and  Dryer  and 
Glassman  (1973)  in  both  carbon  monoxide  and  methane  oxidatioiL  Westbrook  et  al. 
(1977)  developed  a  detailed  finite  rate  model  to  analyze  the  experimental  results  of 
Dryer  and  Glassman.  The  resulting  mechanism  consisted  of  19  species  and  56  reactions 
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which  was  validated  for  the  temperature  range  1000-1350  K.  A  subset  of  this  reaction  set 
was  presented  by  Dryer  and  Glassman  (1978)  for  Hj  -CO  oxidation.  Correa  et  al.  (1984) 
presented  a  partial  equilibrium  model  for  a  turbulent  CO-Hj  -Nj  coaxial  jet  reacting 
with  air  at  atmospheric  pressure.  The  model  was  an  extension  of  the  two-scalar  pdf 
approach  of  Janicka  and  Kollmaim  (1979)  to  include  CO  in  the  radical  pool.  White, 
Drummond  and  Kumar  (1987)  used  a  double  flame  sheet  model  for  temperatures  below 
2500  K  in  a  dual  combustor  ramjet  analysis  which  considered  and  CO  in  the  fuel. 
Above  2500  K  a  chemical  equilibrium  calculation  was  performed.  The  solution 
procedure  was  based  on  an  explicit  forward  marching  boundary-layer  approach. 

The  first  phase  of  the  present  effort,  involved  application  of  the  CMINT  computer 
code  (Scientific  Research  Associates  1991)  to  both  the  projectile  forebody  flow  and  the 
projectile  base  flow  analysis  both  with  and  without  injection.  Both  an  algebraic  mixing 
length  and  a  two-equation  k-«  turbulence  model  were  employed  in  the  initial  studies. 
The  Baldwin-Lomax  (1978)  model  as  described  by  Sahu  and  Danberg  (1986)  was 
subsequently  implemented  for  the  projectile  forebody  turbulence  model,  and  a  wake 
mixing  model  due  to  Chow  (1985)  was  used  downstream  of  the  projectile  base. 

Subsequently,  several  combustion  models  which  are  applicable  to  the  projectile 
base  burning  flow  problem  were  developed.  Application  of  these  models  demonstrates 
the  effect  of  base  region  burning  on  the  projectile  base  pressure.  In  addition  to  these 
projectile  flows,  a  validation  calculation  was  performed  for  comparison  with  the 
experimental  data  of  Jarrett  et  al.  (1988)  on  a  supersonic  biumer  (SSB)  using  fuel. 


2.  ANALYSIS 

The  present  combustion  model  development  effort  focused  on  finite  rate  reaction 
models  which  were  general  enough  to  encompass  the  flow  conditions  encountered 
throughout  the  flight  regime  of  current  and  proposed  Army  base  burning  projectiles. 
Since  the  flow  behind  a  projectile  contains  recirculation  zones,  the  reaction  schemes 
considered  must  be  suitable  for  inclusion  in  a  Navier-Stokes  analysis.  An  implicit 
numerical  procedure  is  desirable  because  of  both  the  presence  of  thin  shear  layers  and 
the  probable  stiff  nature  of  the  equations  due  to  the  chemical  source  terms  in  the 
species  conservation  equations. 
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2.1  Governing  Equations.  The  equations  describing  the  viscous,  chemically 
reacting  projectile  base  flow  are  the  ensemble-averaged  Navier-Stokes  equations 
coupled  with  the  species  conservation  and  turbulence  model  equations.  The  mean  flow 
equations  are  obtained  by  using  mass-weighted  (Favre)  averages  of  the  dependent 
variables.  For  the  present  application  these  equations  are  written  in  a  nonorthogonal 
body-fitted,  t^lindrical  coordinate  system.  The  governing  partial  differential  equations 
were  formulated  in  conservation  form  by  application  of  a  Jacobian  transformation  to  the 
equations  in  ^lindrical  coordinates.  An  outline  of  the  transformation  as  well  as  the 
transformed  system  of  equations  is  given  in  Appendix  A.  The  vector  form  of  the 
equations  is  described  below. 

The  continuity  equation  is  written  as 

+  V.  (pXS)  =  0  (1) 


The  momentum  conservation  equation  is 


3(pXJ) 

at 


+  V.  (pun) 


~Vp  +  V*T 


(2) 


where  r  is  the  stress  tensor  (molecular  and  turbulent)  given  by 


’’ij  =  2peff  ©ij  -  -  peff  v-u  fiij 


(3) 


and  the  rate  of  strain  tensor,  ejj  is  given  by 


1  r  aUj  aUj  . 

2  [  axj  axi  J 


(4) 


The  effective  viscosity,  is  the  stun  of  the  molecular  and  turbulent  viscosities 


f^Bff  ^  P  +  Pt  (5) 

The  turbulent  viscosity,  pj,  is  obtained  from  the  turbulence  model. 

The  energy  conservation  equation  is  written  in  terms  of  the  stagnation  enthalpy, 
ho, as 


3  (phg  j 


at 


+  V 


(oOh.)  -  If  - 


V* q  +  (t •  U) 


(6) 
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where  the  last  term  in  Eq.  (6)  is  the  stress  work  and  q  is  the  multicomponent  energy  flux 
vector  consisting  of  the  Fourier  heat  flux  and  interdiffusional  energy  flux  q^, 


9  =  -  *eff  +  q<j  (7) 

where  is  the  effective  thermal  conductivity.  In  the  present  analysis,  k  is  obtained 
assuming  constant  molecular  and  turbulent  Prandtl  numbers,  Pr  and  Prj,  i.e. 


/iCp  /iTCp 

'  'Ceff  =  -TT  + 

Pr  Pr^ 

The  interdiffusional  energy  flux  is  given  by 
Ns 

=  X  hi(T)  ji 

i=l 


(8) 


(9) 


where  is  defined  in  Eq.  (12)  and  hi(T),  the  enthalpy  of  species  i  per  unit  mass,  is 


r’’ 

hi(T)  =  hfi  +  J^^Cpi(T')dT' 
The  species  conservation  equations  are  expressed  as 


(10) 


+  V.  (pOYi)  .-v.Ji  +  mi 


(11) 


where  Yj  is  the  mass  fraction  of  species  i,  m|  is  rate  of  production  of  species  i  due  to 
chemical  reaction,  and  j;  is  the  diffusional  mass  flux  of  species  i.  Assuming  that  the 
diffusion  of  mass  is  governed  by  Pick’s  law,  j|  is  given  by 


ji  =  -  pD  VYi  (12) 

where  D  is  the  diffusion  coefficient  (independent  of  species  i)  which  is  obtained  by 
assuming  constant  molecular  and  turbulent  Schmidt  numbers.  Sc  and  SO]-,  i.e. 


ii 

pD  =  —  +  - 

Sc  SCcji 

Finally,  for  a  mixture  of  perfect  gases  the  equation  of  state  is 


(13) 
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p  =  pRT 


(14) 


Ns 

I 

i=l 


Yi 

Wi 


where  is  the  universal  gas  constant,  W;  is  the  molecular  weight  of  species  i,  and  is 

the  total  number  of  species  in  the  system.  The  caloric  equation  of  state  relates  the 
temperature  and  the  static  enthalpy  as 

I 

Ns 

h  -  X  Yi  hi(T)  (15) 

i*l 

This  relation  is  evaluated  using  the  JANNAF  database  of  polynomial  curve  fit 
coefficients  for  Cp^  and  hj  as  functions  of  T  v^ch  are  available  from  NASA  Lewis 
Research  Center  (Gordon  and  McBride  1976). 


22  General  Chemistry  Model.  The  typical  solid  propellant  used  in  base 
burning  projectiles  is  a  fuel  ridi  mixture  which  yields  combustion  products  consisting 
primarily  of  H,,  CO,  HQ,  CO,,  H,0  and  N,.  The  mole  fractions  of  these  constituents 
m  chemical  equilibrium  sum  to  0.997,  hence  there  is  little  error  in  ignoring  the  remaining 
trace  species.  Since  the  available  energy  in  the  HQ  is  relatively  small  compared  to  that 
of  H,  and  CO,  the  HQ  has  been  replaced  a  combination  of  CO,  CO,  and  N, .  Both 
the  heat  of  combustion  and  the  molecular  weight  of  the  equivalent  mixture  were 
matched  to  those  of  the  original  equilibrium  combustion  products.  The  composition  of 
this  equivalent  mixture  is  given  in  Table  L 

As  the  base  injectant  gas  mixes  with  the  free  stream  air,  further  reaction  takes 
place  in  the  region  near  the  projectile  base.  Exactfy  viiere  the  combustion  occurs  is  a 
function  of  the  injectant  gas  temperature,  the  mass  and  momentum  flow  rates,  the 
degree  of  turbulent  mixing,  the  effect  of  turbulent  fluctuations  and  the  rates  of  the 
important  chemical  reactions.  In  the  absence  of  turbulence,  the  reaction  rates  are  fairly 
well  known  for  the  H,  -CO  ^tem;  however,  there  are  still  some  uncertainties  which 
must  be  recognized  in  evaluating  the  results.  The  sensitivity  of  the  reaction  rate 
coefficients  to  pressure  (e.g.,  Gardiner  1984)  does  not  iq>pear  to  be  significant  for  the 
range  of  pressures  encountered  in  the  projectile  base  flow  problem. 

The  injectant  gas  temperature  can  be  calculated  if  it  is  assumed  that  the  solid 
propellant  combustion  products  are  in  chemical  equilibrium  upon  exiting  the  projectile 
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base.  However,  available  thermocouple  measurements  of  the  temperature  in  the 
combustion  chamber  of  a  typical  Army  projectile  (Kayser,  Kuzan  and  Vasquez  1987) 
indicate  a  temperature  (T^qj  «  1500  K)  about  500  K  lower  than  the  predicted  adiabatic 
equilibrium  flame  temperature  obtained  using  the  NASA  Lewis  CET86  code  (Gordon 
and  McBride  1976).  Therefore,  in  the  present  analysis  the  composition  of  the  injectant 
gas  was  determined  by  assuming  that  the  products  were  in  equilibrium  at  the 
experimentally  observed  temperature.  It  should  be  noted  that  the  differences  between 
the  compositions  in  these  two  situations  were  not  large.  Hence,  an  important  factor 
regarding  the  injection  conditions  is  simply  the  correct  specification  of  the  injectant  gas 
temperature  because  of  its  direct  influence  on  base  pressure. 

In  general  reactions  are  of  the  form 


I 


I 

i=l 


v' 

i/r 


^f,r 


i 


I 


i=l 


i/r 


(16) 


where  and  v",  are  the  stoichiometric  coefficients  appearing  on  the  left  and  right  of 
the  reaction  r  and  and  Iq,^  are  the  forward  and  backward  rate  constants  and  [Xj]  is  the 
molar  concentration  of  any  species  Xj.  It  can  be  shown  (e.g.,  Vincenti  and  Kruger  1965) 
that  the  rate  of  production  of  species  i  in  reaction  r  is  given  by 

i 

n.r  1  kf.r  H  W 

''  S=1 


where  the  forward  and  backward  reaction  rates  are  related  by  the  equilibrium  constant, 
Kc^viz., 


Kc,r  =  ^  <«) 

Reaction  rates  for  the  models  considered  herein  are  expressed  in  the  Arrhenius  form, 
i.e.. 
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(19) 


l=f,r  ■=  Ar  T*=r  exp(-  ^  ] 

where  R„  is  the  universal  gas  constant,  and  Ej.  is  the  activation  energy  of  reaction  r. 

The  equilibrium  constant  is  actually  a  function  of  temperature  and  is  given  by  the 
relationship  (Vincenti  and  Kruger  1965), 


(20) 


where  yil  is  given  by 

T  T  C*N 

2J(T)  =  [  Sp  dT  +  fil  -  T  I  [  dT  +  Si  I  -  Ru  in  p. 

•■To  L  J  (21) 

AO  AO 

where  hj  and  Sj  are  the  enthalpy  and  entropy  per  mole  at  the  reference  conditions 
and  pg.  The  quantity  is  the  chemical  potential  (Ciibbs  free  energy  per  mole)  for  a 
pure  species  at  unit  pressure.  Eqs.  (18)  -  (21)  enable  the  calculation  of  both  the 
forward  and  backward  rate  constants,  provided  that  the  constants  k^  A,,  b,  and  are 

known  for  each  reaction,  and  Eq.  (17)  represents  the  rate  source  term  for  each  species. 
Since  the  species  equations  are  written  with  the  mass  fractions  as  the  dependent 
variables,  the  molar  concentrations  are  related  to  the  mass  fractions  by 

pYi 

M  =  (22) 

Thus,  the  source  terms  due  to  both  forward  and  backward  reactions  can  be  expressed  in 
terms  of  the  dependent  variables:  the  three  velocity  components,  the  density,  the 
stagnation  enthalpy,  and  the  mass  fraction;  and  can  be  appropriately  linearized  for 
implicit  treatment.  The  influence  of  both  concentration  fluctuations  on  the  chemical 
production  rates  and  temperature  fluctuations  on  the  Arrhenius  rate  expressions  has 
been  neglected  in  this  study. 

The  ignition  of  the  injectant  gas  should  be  quite  rapid  under  projectile  base 
injection  conditions  at  projectile  launch,  since  the  temperature  is  high  and  the  pressure  is 
near  one-half  atmosphere;  therefore  reaction  mechanisms  that  properly  model  the 
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chemistry  of  the  ignition  process  may  not  required.  The  sensitivity  of  the  analysis  to 
the  initial  base  pressure  must  be  investigated  for  extremes  in  the  initial  flight  conditions. 
The  first  mechanism  considered  for  CO  oxidation  is  due  to  Westbrook  et  ai.  (Dryer  and 
Classman  1978);  however,  the  rates  presented  may  not  be  valid  for  temperatures  in 
excess  of  1350  K.  Therefore,  rates  have  been  selected  from  the  work  edited  by  Gardiner 
(1984)  wherever  possible.  The  modified  model  reactions  and  the  Arrhenius  constants 
for  the  forward  reactions  are  shown  in  Table  n  for  reaction  set  A,  This  set  consists  of  23 
reactions  involving  nine  primary  species  ,  O, ,  OH,  Hj  O,  O,  H,  CO,  COj ,  and  Nj  in 
addition  to  two  secondary  species  HO^  and  H^O, ;  this  set  should  permit  the 
calculation  of  ignition  and  the  evaluation  of  its  importance  for  projectile  base 
combustion  flows. 

Since  the  initial  results  for  the  base  combustion  problem  with  set  A  indicated  very 
rapid  ignition  due  to  the  high  initial  injection  temperature,  a  reduced  set  of  species  and 
reactions  could  be  considered.  The  effect  of  reaction  set  A  on  the  projectile  base  drag 
may  be  neglected  compared  to  the  reduced  reactions  sets  considered  under  this  study 
since  a  maximum  change  of  0.002  in  Cdb  values  was  observed.  The  reactions  comprising 
set  B  are  the  first  12  reactions  involving  the  nine  primary  species.  The  first  eight 
reactions  in  Table  11  are  the  same  as  the  Spiegler,  Wolfshtein  and  Manheimer>Tiinnat 
(1976)  reaction  set  Reaction  number  nine  was  also  used  by  Eklund,  Drummond  and 
Hassan  (1990),  but  the  rate  proposed  by  Westbrook  et  al.  (1977)  is  much  larger  than  that 
used  by  Eklund,  Drummond  and  Hassan  (1990).  In  addition,  the  work  of  Gardiner 
(1984)  quotes  a  rate  between  the  latter  two,  but  with  unspecified  products  of  reaction. 
The  importance  of  this  reaction  has  been  assessed  for  the  present  applications. 

23  Global  Hydrogen  -  Air  Combustion  Model.  In  the  case  of  pure  H, 
combustion  in  air,  some  simplification  of  the  finite  rate  chemistry  model  is  possible 
under  certain  circumstances.  When  the  prediction  of  ignition  is  not  a  critical  factor,  say 
due  to  high  initial  temperatures,  then  the  Rogers  and  Chinitz  (1983)  two-step  global 
reaction  model  is  very  attractive.  This  model  requires  only  five  species  including  N, ; 
hence,  it  is  more  efficient  than  the  detailed  reaction  mechanisms.  The  model  is  strictly 
applicable  only  for  combustion  at  one  atmosphere  pressure  and  the  effect  of 
stoichiometry  on  the  global  reaction  rates  is  included.  This  model  was  selected  for  the 
base  combustion  calculation  in  order  to  demonstrate  a  simpler  reaction  set  than  that  for 
H,  -CO  combustiotL  The  applicability  of  the  model  is  assessed  by  comparison  with  other 
base  combustion  model  results. 
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The  reactions  for  the  Rogers  and  Chinitz  model  are: 


+  Oj  ~  20H  (23) 

and 

20H  +  H,  zzr  2H3O  (24) 

where  in  this  form  the  Arrhenius  equation  is 

kf^  «  Ai(<^)T‘'i  (25) 

where  ^  is  the  equivalence  ratio  for  the  overall  reaction  process.  For  the  first  reaction 
(.^)  =  (8.917«  +  3.433/^  -  28.950)  x  10^^ 

E^  •=  4865  cal/nole  (26) 

bj  =  -10 

The  dimensions  of  kf^  are  cm’ /mole-sec.  For  the  second  reaction 

Aj(«)  =  (2.000  +  1.333/«^  -  0.833^)  x  10*< 

E,  =  42,500  cal/mole  (27) 

bj  =  -13 

The  dimensions  of  kf  ^e  cm* /mole  ^ -sec  The  equivalence  ratio  is  defined  as: 


(F/0) 

{k  ss  ■ 

^  (F/0) St 

where  (F/0)  is  the  fuel  to  oxidizer  ratio  by  mass. 


(28) 
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For  the  reaction  set  of  Eq.  (23)  and  Eq.  (24) 


=  5 


1/8 


(29) 


2.4  Turbulence  Models.  Both  an  algebraic  mixmg  length  model  and  a 
two-equation  k-e  model  due  to  Jones  and  Launder  (1972)  were  previously  included  in 
the  CMINT  code.  In  addition,  the  Baldwin-Lomax  (1978)  model  as  described  by  Sahu 
and  Danberg  (1986)  was  implemented  for  the  projectile  forebody  turbulence  model,  and 
a  wake  mixing  model  due  to  Chow  (1985)  was  used  downstream  of  the  projectile  base. 
The  k-e  model  was  evaluated  in  some  preliminary  base  flow  calculations. 

The  approach  taken  in  these  models  assumes  an  isotropic  turbulent  viscosity,  tij, 
and  relates  the  Reynolds  stress  tensor  to  the  mean  flow  gradients,  viz. 

-  p  u'u'  =  pt  [2«ij  “  I  ^ij] 

where  ejj  is  given  in  Eq.  (4). 

2.4.1  Algebraic  Mixing  Length  Model.  In  the  mixing  length  model  the  turbulent 
viscosity  is  determined  from 

Pt  =  (31) 

and  the  mixing  length  is  obtained  from  the  Van  Driest  formulation  with  a  free  stream 
length  scale  Jt„ 


=  0.09  S 


(32) 


where  £  is  the  local  boimdaiy  layer  thickness.  The  mixing  length  is  given 
1  =  tanh 


(33) 


where  is  the  distance  normal  to  the  wall  and  /c  is  the  von  Karman  constant,  k  -  0.4. 


11 


The  van  Driest  damping  factor  Djj  is 


Dx  =  1  -  e 

where  the  wall  coordinate  y  is 


-yV26 


=  —  Yn 


and  the  friction  velocity  is  obtained  from  the  wall  shear,  ,, 

rTwi  h 


=  K] 


(34) 


(35) 


(36) 


While  this  model  gives  acceptable  results  for  turbulent  viscosity  on  the  projectile 
forebody,  it  does  require  determination  of  the  boimdaiy  layer  thickness.  Thus,  this 
model  is  not  directly  applicable  to  the  base  region  turbulence. 


2.42  Baldwin-Lomax  Model.  The  present  description  of  the  so-called 
Baldwin-Lomax  (1978)  model  was  summarized  Ity  Sahu  and  Danberg  (1986).  First,  an 
"inner"  layer  turbulent  viscosity  is  defined  by 


(/^T)inner  =  (37) 

where  the  Van-Driest  mixing  length  is  given  by 

^  =  *ynDx  (38) 

and  the  damping  factor  is  defined  in  Eq.  (34).  Also,  |  (•>  |  is  the  absolute  value  of  the 
vorticity.  The  "outer”  layer  turbulent  viscosity  is  defined  by 

(^t)  outer  “  ^  ^cp  P  ^wake  ^kleb(y)  (38) 

where 

^wake  “  ®^”(yinax^inax»  ^kymax'^*diff/^niax) 

Fmax  “  max[F(y)]  -  F(yjnax)  (40) 

F(y)  =  ynl"|Dx 

8n<l  is  the  difference  between  the  maximum  and  minimum  velocities  in  a  shear 
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layer.  For  a  boundary  layer,  the  above  minimum  velocity  is  zero.  Also,  the  Van  Driest 
damping  factor  is  neglected  in  free  shear  layers  and  wakes.  The  constants  used  in  the 
model  are  a:  =  0.4,  C^,  =  1.6,  Cyeb  =  0.3,  =  0.25  and  K  =  0.0168. 

2.43  Jones-Launder  k-t  Model.  The  low  Reynolds  number  k-€  model 
of  Jones  and  Laimder  (1972)  does  not  require  the  specification  of  a  length  scale  or 
boundary  layer  thickness.  One  disadvantage  of  the  model  is  the  requirement  for  fine 
near  wall  resolution  to  resolve  large  gradients  in  the  turbulence  kinetic  energy  (k).  Also, 
the  equations  contain  ad  hoc  low  turbulence  Reynolds  number  correction  terms.  With 
the  k-€  model  the  turbulent  viscosity  is  obtained  from  the  Prandtl-Kolmogorov  relation, 

pk^ 

The  empirical  constants  o^.,a^  and  C,  required  in  the  k  and  e  transport  equations 
(Appendix  A)  are  taken  from  Jones  and  Launder  (1972)  and  the  constant  from 


Launder  and  Spalding  (1974),  i.e., 

Cj  =  1.44  (42a) 

Cj  =  1.92  [l.O  -  0.3  exp[-RT*]]  (42b) 

=  0.09  exp[-  ^  ]  (420) 

=  1.0  (42d) 

=  1.3  (42e) 


2.4.4  Eggers  Turbulence  Model.  The  Eggers  (1971)  algebraic  mixing  length  model 
was  developed  for  the  nonreacting  mixing  of  coaxial  hydrogen-air  jets.  For  the  reacting 
coaxial  jet  experiment  of  Jarrett  et  al.  (1988),  this  model  was  modified  by  Eklund, 
Drummond  and  Hassan  (1990)  to  use  the  diatomic  hydrogen  profile  to  characterize  the 
shear  layer.  The  turbulent  viscosity  was  defined  as 


fiiji  —  Cq  p  Ucl  (43) 

where  Q  is  a  constant  (0.032),  Ud  is  the  streamwise  velocity  on  the  jet  centerline,  and  Rj 
is  the  width  of  the  mixing  layer.  The  width  is  defined  as  the  radial  distance  between  the 
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points  in  the  profile  where  the  Hj  mass  fractions  are  (H^ )  and  Y,  (Hj ), 

=  YaCHa)  +  0.95  [Yo(Hj)  - 

^7(^2)  =  ^3(^2)  +  0.5  [Yo(Hj)  -  YaCH,)]  (44) 

where  Y3(H2  )  is  the  H,  mass  fraction  in  the  external  outer  jet  flow  and 
Yo(H,)  is  the  mass  fraction  of  H,  on  the  jet  centerline. 

I 

25  Solution  Technique.  Solutions  of  the  above  equations  were  computed  using  a 
reacting  flow  version  of  SRA’s  Navier-Stokes  code,  CMINT.  Centered  spatial 
differences  were  used  with  adjustable  artificial  dissipation.  The  equations  were  solved 
using  a  linearized  block  implicit  (LBI)  algorithm  and  an  ADI  approximate  factorization. 
A  spatially  varying  time  step  was  used  to  accelerate  convergence  to  a  steady  solution. 

For  the  reacting  flow  solutions  the  time  step  for  the  coupled  species  equations  was 
further  conditioned  using  a  time  step  scaling  based  on  the  chemical  production  source 
terms  in  the  species  equations.  The  sharp  comer  at  the  projectile  base  was  treated  as  a 
grid  cut-out  region  using  a  single  non-rectangular  computational  domaiiL  A  more 
complete  description  of  the  solution  technique  is  given  by  Briley  and  McDonald  (1977, 
1980). 

The  approach  used  by  Eklund,  Drummond  and  Hassan  (1990)  was  to  treat  the 
chemical  source  terms  implicitly  on  a  pointwise  basis.  Since  an  explicit  solution 
procedure  was  used  1^  Eklund,  Drummond  and  Hassan  (1990),  this  amounts  to  rescaling 
the  time  step  in  each  individual  species  equation,  and  therefore  allows  each  equation  to 
relax  at  its  own  time  scale.  In  this  approach  the  species  equations  are  still  solved  in  an 
uncoupled  maimer  at  each  time  step. 

The  present  fully  implicit  approach  automatically  includes  the  pointwise  implicit 
coupling  of  the  chemical  source  terms,  and  the  coupling  among  the  various  species 
equations  being  solved.  The  CMINT  code  allows  for  the  user  specified  coupling  of  the 
species  equations,  and  coupling  to  the  Navier-Stokes  equations  is  optional  as  well. 
Therefore,  if  certain  species  do  not  contribute  significantly  to  the  eneigy  balance,  those 
equations  could  be  solved  decoupled  from  the  other  species  and  flow  equations  at  each 
time  step.  This  approach  can  save  computer  time  per  time  step  although  the 
convergence  rate  may  be  adversely  affected. 

The  coupling  of  the  species  equations  involved  in  energetic  reactions  is  important 
and  improves  overall  convergence.  However,  the  possibility  exists  that  the  species 
equation  source  terms  will  cause  ill-conditioned  matrices  due  to  large  off-diagonal 
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elements  in  the  block  matrix  at  a  particular  grid  point.  Obviously,  this  problem  is  not 
encountered  when  the  species  equations  are  solved  one  at  a  time  (the  decoupled 
approach).  In  order  to  solve  this  difficulty  with  ill-conditioned  block  matrices,  a  time 
step  scaling  was  devised  for  the  coupled  species  equations.  This  scaling  factor  was 
applied  in  addition  to  the  spatial  time  step  conditioning,  which  is  applied  to  all  equations. 

The  system  of  P.D.E.’s  may  be  written  as 


aH(^) 

at 


=  D  <#>  +  S(<^) 


(45) 


where  D  is  the  nonlinear  spatial  difference  operator  matrix  and  S(^)  is  the  source  term  vector. 
The  linearized  difference  equations  are  written  as 


A  ^  +  1^^  +  b(^^) 


(46) 


where 


-n+1  -n 

as  ^ 


(47) 


and 


a-22 

a^ 


aD 

£  =  — 

a^ 


(48) 


The  species  equation  source  terms  would  appear  as  S(0)  terms  and  these  terms  can 
cause  ill-conditioned  matrices  due  to  the  form  of  the  chemical  production  terms  in 
Eq.  (17). 

The  conditioning  factor  for  the  coupled  species  equations  is  chosen  to  insure 
diagonal  dominance  of  the  block  matrix.  First,  the  maximum  absolute  value  of  the 

operator  a  S/a  off-diagonal  elements  is  determined  at  each  grid  point  for  each  of  the 
coupled  species  equations.  Then  the  time  derivative  matrix  A  is  replaced  by 


FfAij 


(49) 
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for  the  coupled  species  equation.  The  factor  Fj  is  given  by 


Fi  =  max 


dSj 
“ax  — 


1.0, 


Aii/At 


(50) 


where  the  constant  Q  is  an  arbitrary  parameter  with  values  of  0.1  and  1.0  used  in  the 
present  applications.  This  amounts  to  rescaling  the  time  step  operators  for  p  and  Yj  in 

I 

each  of  the  coupled  species  equations.  Note  that  this  scaling  is  neither  defined  nor  used 
for  a  decoupled  species  equatioiL 


2.6  Two-Phase  Flow  Analysis.  Two-phase  flow  effects  may  be  present  in  certain 
projectile  base  bum  applications  depending  primarify  on  the  propellant  formulation, 
particle  size  distribution,  and  burning  rate.  The  M864  base  bum  projectile  uses  an 
ammonium  perchlorate  (AP)  oxidizer  based  fuel,  and  generally  is  e}q>ected  to  generate 
mostly  gaseous  combustion  products  upon  exiting  the  projectile  base.  Some  small  AP 
particulates  (^5  pm  or  less)  may  remain  in  the  injected  gas.  Therefore,  a  sophisticated 
two-phase  flow  analysis  for  projectile  base  flow  q)plications  is  probabfy  not  required. 
An  existing  two-phase  flow  code  (Sabnis,  Gibeling  and  McDonald  1987)  vm  adapted  to 
the  projectile  application  and  tested  on  a  representative  base  combustion  problem. 

C^omputational  techniques  used  in  simulation  of  two-phase  flows  can  be  broadly 
categorized  into  two  approaches,  viz.  the  Euleiian-Eulerian  analysis  and  the 
Eulerian-Lagrangian  analysis.  Both  techniques  involve  computing  the  continuous  phase 
using  an  Eulerian  analysis.  The  influence  of  the  discrete  phase  (either  solid  particles  or 
liquid  droplets)  on  the  continuous  phase  is  accounted  for  1^  inclusion  of  inter-phase 
coupling  terms  in  the  Eulerian  equations,  which  in  the  absence  of  these  terms  would  be 
the  usual  Navier-Stokes  equations.  The  discrete  phase,  on  the  other  hand,  may  be 
treated  with  either  a  continuum  model  or  a  discrete  model.  The  Eulerian-Eulerian 
technique  uses  a  continuum  model  for  the  discrete  phase  and  is  commonly  termed  the 
two-fluid  model  This  approach  models  a  dense  granular  bed  very  conveniently  and  this 
undoubtedly  accounts  for  its  popularity  in  modeling  gun  interior  ballistics  where  large 
particle  loading  ratios  occur  over  most  of  the  cyde  (e.g.,  Gough  1977  and  Gibeling, 
McDonald  and  Banks  1983).  The  Eulerian-Lagrangian  approach  employs  a  Lagrangian 
description  to  analyze  the  motion  of  the  discrete  phase,  using  con^)utational  "particles" 
to  represent  a  collection  of  physical  particles.  Newton’s  law  of  motion  is  employed  to 
simulate  the  particle  motion  under  the  influence  of  the  local  environment  produced  by 
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the  continuous  phase.  The  discrete  phase  attributes  (such  as  the  particle  position  and 
velocity  vectors,  size,  temperature,  etc.)  are  updated  along  the  trajectories. 

In  simulation  of  flows  containing  burning  particles  or  evaporating  droplets,  it 
becomes  necessary  to  account  for  the  fact  that  the  discrete  phase  is  not  mono-dispersed. 
To  accomplish  this  in  the  Eulerian-Eulerian  methodology,  the  two-fluid  model  can  be 
generalized  into  a  multi-fluid  model.  However,  the  CPU  time  requirements  increase 
rapidly  with  increasing  number  of  particle  size  classes,  since  an  extra  "fluid"  has  to  be 
added  for  every  particle  size  class,  therely  increasing  the  number  of  partial  differential 
equations.  The  Eulerian-Lagrangian  analysis,  on  the  other  hand,  treats  the  particle  size 
as  one  of  the  attributes  assigned  to  computational  particles  and  hence  has  no  trouble 
simulating  flows  which  involve  changing  particle  size.  Since  this  approach  involves 
integration  of  ODE’s  for  the  particulate  phase,  it  is  numerically  efficient.  Furthermore, 
the  deterministic  nature  of  the  particle  dynamics  facilitates  the  incorporation  of  models 
for  turbulent  dispersion,  agglomeration,  collision,  etc. 

In  Eulerian-Lagrangian  algorithms,  the  inter-phase  coupling  terms  for  the 
continuous  phase  equations  can  be  computed  using  a  particle  trajectory  approach  or  a 
particle  distribution  approach.  In  the  particle  trajectory  approach,  the  coupling  terms 
are  computed  from  the  knowledge  of  the  trajectories  for  representative  particles  and 
their  attributes  at  the  intersection  of  the  trajectories  with  the  Eulerian  cell  boundaries. 

In  the  particle  distribution  approach,  the  coupling  terms  are  computed  ffom  the 
instantaneous  distribution  of  the  particles  in  the  computational  domain.  The  trajectory 
approach  has  been  employed,  for  example,  by  Crowe,  Sharma  and  Stock  (1977),  and 
Gosman  and  loannides  (1983),  while  the  particle  distribution  approach  has  been  utilized 
by  Dukowicz  (1980)  and  Sabnis,  Gibeling  and  McDonald  (1987). 

In  the  algorithms  based  on  the  trajectory  approach,  the  integration  of  the 
Lagrangian  equation  of  motion  for  representative  particles  is  carried  out  starting  from 
the  injection  location  until  the  particle  leaves  the  computational  domain  or  imtil  its  size 
becomes  negligible.  During  this  interaction,  the  continuous  phase  flow  Held  is  held 
frozerL  The  inter-phase  coupling  terms  for  the  continuous  phase  conservation  equations 
are  computed  for  every  Eulerian  cell  from  the  net  influx  of  the  appropriate  conserved 
variable  into  the  Eulerian  cell,  due  to  all  trajectories  intersecting  the  particular  Eulerian 
cell.  The  coupling  terms  thus  computed  are  used  to  calculate  the  continuous  phase  flow 
ff  eld  which  can  then  be  used  to  re-evaluate  the  trajectories  and  the  source  terms.  This 
iterative  process  is  continued  until  the  desired  level  of  convergence  is  achieved.  These 
algorithms  are  thus  inherently  unsuitable  for  transient  calculations  and,  further,  the 
global  iteration  procedure  used  can  require  substantial  computer  time. 
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In  the  particle  distribution  approach,  such  as  that  used  by  Dukowicz  (1980)  and 
Sabnis,  Gibeling  and  McDonald  (1987),  the  source  terms  are  computed  from  the 
instantaneous  interaction  between  the  continuous  phase  and  all  the  particles  in  the 
particular  Eulerian  cell.  Thus,  the  source  term  for  the  continuous  phase  continuity 
equation,  for  example,  is  given  by  the  sum  of  the  mass  transfer  rates  for  all  the  particles 
in  the  cell.  The  calculation  procedure  consists  of  updating  the  particle  distribution 
through  one  time  step  followed  by  updating  the  continuous  phase  flow  field  through  one 
time  step.  In  general,  it  is  not  necessary  that  the  time  step  used  to  integrate  the  particle 
motion  and  that  used  in  updating  the  continuous  phase  flow  field  be  equal.  However,  by 
making  the  two  time  steps  equal,  the  particle  distribution  algorithms  can  be  used  for 
simulation  of  transient  phenomena.  If  only  a  steady-state  solution  is  desired,  then  the 
two  time  steps  can  be  made  unequal  and  matrix  preconditioning  techniques  can  be  used 
for  convergence  acceleration  of  the  continuous  phase  solution. 

The  present  analysis  is  based  on  the  CELMINT  (Combined  Eulerian  Lagrangian 
Multidimensional  Implicit  Havier-Stokes  Jime-dependent)  code  developed  by  Sabnis, 
Gibeling  and  McDonald  (1987).  In  this  algorithm,  the  ensemble-averaged  Navier-Stokes 
equations  (including  the  inter-phase  coupling  terms)  are  solved  for  the  continuous  phase. 
A  particle  distribution  model  is  used  in  the  Lagran^an  treatment  of  the  particulate 
phase.  The  key  feature  of  the  particle  transport  model  in  CELMINT  is  that  it  integrates 
the  Lagrangian  equations  of  motion  for  a  particle  in  conq>utational  space  rather  than 
physical  space.  This  simplifies  the  computation  of  the  interphase  coupling  terms, 
because  the  search  for  the  mesh  cell  location  of  a  particle  becomes  trivial 

The  CELMINT  code  has  been  validated  previously  (cf.  Sabnis  et  al.  1988)  using 
the  experimental  data  reported  by  Milojevic,  Bomer  and  Durst  (1986)  for  two-phase 
shear-layer  flow  without  inter-phase  mass  transfer.  More  recently  (cf.  Sabnis  and  de 
Jong  1990),  this  Eulerian-Lagrangian  analysis  was  utilized  to  simulate  the  two-phase  flow 
in  an  evaporating  spray  and  the  calculated  results  were  compared  with  the  ejq>erimental 
data  of  Solomon  et  al.  (1984).  The  equations  to  be  solved  for  the  continuous  phase  are 
the  mass,  momentum,  and  energy  conservation  equations  including  the  appropriate 
source  terms  to  account  for  the  influence  of  the  particulate  phase  on  the  continuous 
phase.  The  form  of  these  terms  for  a  single  species  particle  is  given  in  Sabnis,  Gibeling 
and  McDonald  (1987)  and  Sabnis  and  de  Jong  (1990). 

Under  the  present  effort  the  Lagrangian  module  was  modified  for  application  to 
projectile  base  combustion  with  particles  using  the  CMINT  code.  This  module  could  be 
implemented  in  other  Navier-Stokes  codes  if  desired.  For  the  base  flow  problem  a 
boundary  definition  routine  is  required  to  permit  calculation  of  particle  motion  with 
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realistic  boundary  interactions.  A  restitution  coefficient  model  is  used  for  particle-wall 
collisions  (de  Jong,  Sabnis  and  McConnaughey  1989).  The  present  application  has  been 
tailored  for  the  analysis  of  ammonium  perchlorate  ( AP)  vaporization,  since  these  are  the 
most  likely  particles  to  be  emitted  from  the  projectile  base  burning  propellant.  The 
equilibrium  products  for  the  self-deflagration  of  AP  are  Oj ,  N, ,  O  and  HQ;  it  is 
reasonable  to  replace  the  HCl  with  an  equivalent  amount  of  CX)  and  N, .  Therefore,  the 
particulate  AP  is  assiuned  to  consist  of  the  following  species  in  the  present  analysis. 


Species 

Mass  Fraction  (fj! 

Molecular  Weight 

0. 

0J68 

31.999 

N, 

0354 

28.013 

HjO 

0353 

18.015 

CO 

0.025 

28.01 

Mixtme: 

1.000 

2538 

A  vaporization  model  based  on  a  Sherwood  number  analysis  for  an  isolated 
spherical  particle  has  been  incorporated  into  the  Lagrangian  calculation  procedure.  A 
linear  regression  burning  rate  has  been  used  in  the  analysis,  and  the  burning  rate  has 
been  obtained  from  the  available  AP  strand  burning  e}q>erimental  data  for  the  M864 
propellant.  The  resulting  rate  of  gas  mass  production  of  species  i  due  to  particle  biuning 
under  these  assumptions  may  be  written  as, 

nii  =  fi  (51) 

where  fj  is  the  AP  species  mass  fraction  from  the  above  table.  The  particle  vaporization 
rate  is  assumed  to  enhanced  by  gas  convection  around  the  particle,  and  is  given  by, 

iv  =  -0.5  Sh  [411  Rp  pp]  Rp^t  (52) 

The  rate  of  change  of  particle  radius,  Rp  (,  is  a  negative  constant  for  linear  regression, 
and  the  Sherwood  number,  Sh,  is  the  mass  transfer  analogy  of  the  Nusselt  number  for 
heat  transfer.  The  Sherwood  number  is  assumed  to  be  the  same  as  the  Nusselt  number 
for  an  isolated  spherical  particle,  which  is  based  on  the  relative  velocity  between  the  gas 
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and  particle,  i.e. 


Sh  =  2  +  0.53  (Rep)®*‘  for  RCp  <  21B.92 
Sh  =  0.37  (Rep)“*«  for  Rep  ^  278.92 


(53) 


where  the  particle  Reynold’s  number  is  defined  as. 


I 


Re, 


p  2Rp  |U  -  Upl 


(54) 


Complete  details  of  the  combined  Eulerian-Lagrangian  procedure  using  the 
CMINT  code  are  given  in  Sabnis,  Gibeling  and  McDonald  (1987),  Sabnis  et  al.  (1988), 
and  Sabnis  and  de  Jong  (1990)  and  are  not  repeated  here.  A  sample  calculation  has  been 
performed  and  is  discussed  in  the  section  on  Base  Flow  Applications. 


3.  REACTING  FLOW  VAUDATION  CASE 

A  supersonic  flow  coaxial  burner  (SSB)  studied  ejqierimentally  by  Jarrett  et  aL 
(1988)  has  been  selected  as  a  validation  case  for  the  present  analysis.  While  this  case 
only  considers  H,  combustion,  it  is  well  documented  and  a  digital  version  of  the  data  is 
available  from  Jarrett  et  aL  (1988).  Also,  this  case  has  been  analyzed  numerical^  by 
Jarrett  et  aL  (1988)  and  Eklund,  Drummond  ^d  Hassan  (1990).  A  schematic  of  the  SSB 
apparatus  and  computational  domain  is  shown  Fig.  1.  The  SSB  consists  of  an  inner 
hydrogen  jet  exiting  at  M  « i.O  with  a  coaxial  vitiated  air  jet  exiting  at  M = 2.0.  The 
inflow  boundary  conditions  for  the  calculations,  based  on  the  ideal  burner  exit  conditions 
obtained  from  Jarrett  et  aL  (1988),  are  given  in  Table  IIL  The  SSB  nozzle  walls  are 
conical  with  a  half-angle  of  43  degrees.  The  fuel  injector  is  a  cone-blinder  geometry  as 
shown  in  Fig.  1,  and  a  shock  wave  emanating  from  the  cone-cylinder  juncture  leads  to 
some  uncertainty  in  the  jet  conditions  specified  in  Table  m,  as  noted  Ity  Eklund, 
Drummond  and  Hassan  (1990).  The  Eggers  turbulence  model  as  modified  by  Eklund, 
Drummond  and  Hassan  (1990)  was  used  for  this  case  (see  section  2.4.4). 

Three  different  (Cartesian  grid  systems  have  been  used  on  this  case  to  determine 
the  effect  of  mesh  refinement  on  the  solutiorL  The  first  grid  utilized  101  radial  points 
and  61  axial  points;  the  second  used  101  radial  and  101  axial  points  (Fig.  2);  and  the  third 
used  111  radial  and  121  axial  points.  All  grids  used  nonuniform  distributions  in  both 
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directions.  It  can  be  seen  that  the  lip  between  the  fuel  injector  and  coaxial  air  stream  is 
well  resolved  in  the  second  grid,  while  the  outer  nozzle  lip  has  poorer  resolution.  The 
third  grid  was  constructed  based  on  the  solution  using  the  second  grid  to  better  resolve 
regions  of  steep  gradients. 

Calculations  were  Hrst  made  on  the  three  grids  using  a  nine  reaction  set  consisting 
of  reactions  one  through  nine  from  Table  n.  A  final  calculation  was  made  on  the  third 
grid  deleting  the  ninth  reaction  to  determine  its  importance  in  this  case,  which  was  not 
significant.  All  of  the  calculations  were  started  by  assuming  the  unmixed  jets  extended  to 
the  outflow  boundary  with  a  blending  region  between  the  fuel  and  air  streams.  The 
initial  constant  pressure  throughout  was  set  equal  to  the  pressure  of  the  vitiated  air 
stream.  The  pressures  at  the  hydrogen  exit  and  all  ambient  boundaries  were  modified 
over  100  iterations  (time  steps)  to  achieve  the  values  specified  in  Table  m. 

The  results  of  the  calculation  on  the  second  grid  (101  x  101)  are  shown  in  Figs.  3 
through  6,  and  the  temperature  prediction  on  the  first  grid  (101  x61)  is  shown  in  Fig.  7. 
The  computed  solution  using  the  third  grid  is  very  close  to  that  in  Figs.  3  through  6, 
hence  those  results  are  omitted  here.  The  results  shown  are  for  axial  stations  at  2S.4  mm 
(one  inch)  intervals  starting  at  the  nozzle  exit  The  inflow  axial  velocity  shown  in  Fig.  3a 
indicates  a  significant  difference  in  the  starting  values  used  in  the  CFD  simulations 
versus  the  experimental  results.  This  may  be  caused  experimental  error  due  to  seed 

particle  "lagging"  in  the  high  shear  regions  or  to  distortion  of  the  actual  velocity  profile 
due  to  the  nozzle  and  fuel  injector  configuration.  Also,  as  noted  by  Jarrett  et  aL  (1988) 
the  CARS  and  LDV  measuring  volumes  are  not  small  compared  to  the  fuel  injector 
diameter,  which  will  result  in  flattened  experimental  profiles  in  regions  of  large 
gradients. 

The  velocity  profiles  and  O,  and  N,  number  density  profiles  were  not  significantly 
different  as  a  result  of  the  grid  refinement,  hence  those  figures  for  the  first  grid  have 
been  omitted  here.  The  computed  axial  velocities  are  seen  to  lead  the  experimental 
values  slightly  at  all  measuring  stations,  and  to  a  lesser  extent  in  the  results  of  Jarrett  et 
al.  (1988).  Eklund,  Drummond  and  Hassan  (1990)  did  not  show  velocity  predictions.  It 
should  be  noted  that  the  LDV  turbulence  measurements  Jarrett  et  al.  (1988)  show  large 
anisotropic  turbulent  stresses  which  are  not  modeled  by  simple  algebraic  turbulence 
models  employed  in  the  various  calculations.  Since  much  of  the  initial  flow  field  change 
is  shear  driven,  the  use  of  an  isotropic  model  should  result  in  some  differences  between 
computation  and  experiment 

The  temperature  comparison  with  data  is  shown  in  Figs.  4  and  7,  where  it  is  seen 
that  the  present  results  underpredict  the  core  region  temperature  at  an  axial  location 
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X  =  25.4  mm,  while  over  predicting  the  temperature  somewhat  at  x  =  76.2  mm.  The 
results  on  the  finer  grid  (Fig.  4)  agree  very  well  with  the  data  elsewhere.  The  calculation 
of  Jarrett  et  al.  (1988)  shows  a  similar  discrepanQr  in  temperature  at  x  =  25.4  mm  and  a 
slight  overprediction  at  x  =«  101.6  nun,  and  shows  close  agreement  at  the  other  two 
stations.  Eklund,  Drummond  and  Hassan  (1990)  underpredicted  the  temperature  in  the 
core  at  all  stations  except  x  =  50.8  mm  where  their  results  are  very  close  to  the  data. 

The  data  at  x  =  25.4  mm  possibly  indicates  that  fuel  ignition  has  taken  place  sooner 
than  predicted,  which  is  the  opposite  of  what  is  e]q)ected.  In  a  recent  private 
communication,  Jarrett  (1991)  indicated  the  discovery  of  a  ^tematic  error  in  the  data 
reported  in  Jarrett  et  al.  (1988).  Also,  the  availability  of  more  recent  measurements  on 
the  same  apparatus  Cheng  et  aL  (1991)  was  noted.  The  data  of  Cheng  et  al.  (1991)  is  not 
yet  available  in  digital  form;  however,  the  figures  in  Cheng  et  aL  (1991)  show  that  the 
fuel  has  not  ignited  at  25.4  mm,  and  in  fact  the  present  tenq)erature  prediction  at  that 
location  is  closer  to  this  new  data. 

The  present  predictions  show  somewhat  larger  differences  in  O,  number  density 
(Fig.  5)  than  either  Jarrett  et  al.  (1988)  or  Eklund,  Drummond  and  Hassan  (1990),  while 
the  N,  number  density  is  much  closer  to  the  data.  The  spreading  rate  evident  from  Figs. 
5  and  6  is  somewhat  larger  than  that  obtained  by  either  Jarrett  et  al.  (1988)  or  Eklund, 
Drummond  and  Hassan  (1990).  In  general,  the  level  of  agreement  between  the  present 
predictions  and  experiment  is  quite  good  considering  the  uncertainties  and 
approximations  involved. 

This  validation  case  provides  a  level  of  confidence  in  the  finite  rate  chemistry 
model  implementation  in  the  present  code.  Also,  the  reaction  set  utilized  in  this  case  is  a 
subset  of  the  H,  -CO  reaction  set  used  in  the  base  combustion  calculations,  and  this  case 
implies  a  limited  validation  of  the  reaction  set  and  rate  constants  employed. 


4.  PROJECnUEAPPUCATIONS 

4.1  Boundary  Conditions.  Since  only  supersonic  flow  (M«, »  2)  was  considered  in 
the  present  base  flow  calculations,  the  upstream  boundary  conditions  were  obtained 
from  a  full  projectile  calculation  (Nietubicz  and  Heavey  1990).  Specified  values  for  all 
the  dependent  variables  were  set  on  this  boundary.  For  the  full  projectile  calculations, 
specified  values  were  set  for  the  dependent  variables  on  the  freestream  boundary  ahead 
of  the  projectile.  On  the  outer  radial  boundary  specified  supersonic  conditions  were  set 
from  the  upstream  boundary  to  the  axial  station  of  the  projectile  base,  and  downstreari 
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of  this  station  extrapolation  was  used.  The  outer  boundary  was  located  sufficiently  far 
from  the  projectile  so  that  waves  emanating  from  the  body  pass  through  the  downstream 
boundary.  At  the  projectile  surface  no-slip  conditions,  a  specifled  wall  temperature  (T^ 
=  T„  =  294  K),  zero  normal  pressure  gradient  and  zero  gradient  of  species  mass  fractions 
were  specified.  Along  the  base  injection  region  stagnation  temperature,  axial  mass  flux 
and  species  mass  fractions  were  specified,  while  the  pressure  was  determined  from  the 
normal  momentum  equation  and  the  radial  velocity  component  was  assumed  to  be  zero. 
At  the  downstream  supersonic  outflow  boundary,  first  derivative  extrapolation  was  used. 

4.2  Flat  Base  Projectile  Case.  The  M864  projectile  with  a  flat  nose  and  a  flat  base 
was  considered  to  obtain  a  forebo<fy  flow  field  solution  as  a  starting  condition  for  the 
supersonic  base  flow  computation.  The  projectile  schematic  is  reproduced  in  Fig.  8  from 
Danberg  (1990).  An  algebraic  grid  was  generated  for  this  configuration  with  clustering 
near  the  nose  and  the  projectile  surface  (Fig.  9).  The  resolution  downstream  of  the  base 
was  sacrificed  since  only  the  forebody  solution  was  required  fi'om  this  calculation.  In 
fact,  the  results  shown  were  obtained  by  a<Kuming  an  extended  sting  downstream  of  the 
base.  In  this  case  the  axial  direction  grid  line  emanating  fi’om  the  projectile  base  comer 
was  assumed  to  be  a  solid  surface,  and  points  between  the  centerline  and  this  surface 
were  excluded  from  the  calculation.  The  grid  consisted  of  100  radial  points  and  280  axial 
points  with  28  points  between  the  centerline  and  the  projectile  nose  comer,  65  axial 
points  upstream  of  the  nose,  and  250  axial  points  fiom  the  inflow  boundary  to  the  end  of 
the  projectile. 

The  calculation  was  run  assuming  that  T^  =  T„  =  294  K,  M„  =  2.0,  and  a 
freestream  pressure  of  100,000  Pascals.  A  time  step  conditioning  scheme  (i.e.,  spatially 
varying  time  step)  was  utilized  successfully  to  obtain  good  convergence  throughout  the 
flow  field.  The  siuface  pressure  distribution  for  this  case  is  shown  in  Figure  10,  where 
the  projectile  nose  is  located  at  X/D  =  0.1197.  The  present  results  indicate  a  stagnation 
pressure  on  the  nose  of  Pn/p*  =  5.89,  with  a  rapid  drop  in  pressure  around  the  sharp 
comer.  An  overexpansion  to  a  normalized  pressure  of  034  is  observed,  followed  by  a 
rapid  recover?  iz  the  expected  level.  Comparison  of  the  present  results  with  those  of 
BRL  show  the  two  calculations  to  be  virtually  identical  beginning  slightly  downstream  of 
the  flat  nose  except  for  a  slight  axial  shift,  as  shown  in  Fig.  10. 
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5.0  BASE  FLOW  APPUCATIONS 


5.1  Non-Reacting  Flow  Cases.  Under  the  present  study  the  initial  base  flow 
calculations  were  performed  with  an  algebraically  generated  grid.  A  grid  consisting  of 
130  radial  by  230  axial  mesh  points  was  constructed  for  the  base  flow  calculation.  The 
upstream  boundary  for  this  grid  was  placed  at  X/D  —  435  which  is  about  one  caliber 
upstream  of  the  start  of  the  boattail.  There  were  58  grid  points  between  the  centerline 
and  the  projectile  base  comer.  The  grid  was  reflned  in  the  axial  direction  adjacent  to  the 
base  with  160  points  between  the  base  and  the  downstream  boundary  located 
approximate^  43  calibers  from  the  projectile  base.  Initial  conditions  for  the  base  flow 
calculations  were  obtained  from  the  previously  discussed  forebody  solutioiL  The  base 
flow  calculation  was  started  using  a  simple  mixing  length  model  and  run  for  100  time 
steps  to  allow  the  base  flow  to  develop.  Then  the  calculation  was  restarted  using  the 
Jones  and  Launder  (1972)  low  Reynolds  number  form  of  the  k-c  turbulence  model 

The  calculation  was  run  to  convergence  with  a  reduced  artificial  dissipation 
coefficient  eveiywhere  except  near  the  projectile  base  comer.  A  drag  coefficient 
Cdb  =“  0*  149  was  obtained.  The  base  pressure  distribution  for  this  case,  shown  in  Fig.  11, 
is  relatively  fiat;  however,  the  turbulent  viscosities  predicted  very  near  the  projectile  base 
were  unrealistically  large.  This  may  be  due  to  inadequate  axial  grid  resolution  near  the 
projectile  base,  since  the  Jones-Launder  k-e  formulation  is  sensitive  to  near  wall  grid 
spacing.  Figure  11  also  shows  the  prediction  with  the  BRL  axisymmetric  code  and  a 
mixing  length  turbulence  model  (Baldwin-Lomax/Chow),  and  a  second  k-c  model  result 
using  the  BRL  grid  for  the  base  region  (139  x  166).  This  grid  has  an  axial  spacing  at  the 
base  (axq/D)  of  0.0015  which  is  about  100  times  larger  than  the  spacing  in  the  130  x  230 
grid.  The  latter  grid  shows  some  pressure  oscillations  near  the  comer,  but  otherwise 
only  a  moderate  difference  fi’om  the  more  refined  130  x  230  grid.  As  discussed  in  the 
following  subsection,  52,  when  a  mixing  length  model  was  used  with  the  present  analysis, 
results  which  are  very  similar  to  the  BRL  results  were  obtained. 

A  dome  base  M864  projectile  configuration  was  also  considered  for  the 
nonreacting  cases  considered.  A  base  region  grid  consisting  of  ISO  radial  and  230  axial 
points  was  constmcted  using  the  EAGLE  code  (Thompson  1987).  The  upstream 
boundary  for  this  grid  was  placed  at  X/D  =  435,  and  the  flow  conditions  were  the  same 
as  in  the  fiat  base  case.  The  base  flow  calculation  was  run  with  the  k-c  model  and  a 
reduced  artificial  dissipation  coefficient  except  in  the  immediate  base  regioiL  The  drag 
coefficient  obtained  was  Cdb  =  0'  147.  It  should  be  noted  that  the  results  achieved  with 
the  CMINT  code  showed  sensitivity  to  the  grid  resolution  and  the  second  order  artificial 
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dissipation  coefficient  when  the  k-e  model  was  employed. 

A  further  correction  to  the  k-€  model  was  implemented  to  ascertain  the  influence 
on  base  drag.  The  high  Mach  number  modifications  suggested  by  Childs  and  Caruso 
(1987)  lead  to  a  base  drag  coefficient  of  C^b  =  0.133  for  the  dome  base  case.  The 
streamline  curvature  modification  suggested  by  Childs  and  Caruso  (1989)  was  not 
successfully  employed  in  the  dome  base  case.  However,  the  results  of  Childs  and  Caruso 
indicated  that  the  streamline  curvature  modifications  tended  to  offset  the  high  Mach 
number  corrections  for  the  case  they  considered.  No  additional  consideration  of 
turbulence  model  related  issues  was  carried  out  imder  this  effort 

52  Hot  Injection  and  Reacting  Flow  Cases.  Subsequent  calculations  under  this 
project  were  carried  out  for  a  flat  base  M864  projectile  with  a  three  degree  boattail  angle 
using  a  full  grid  provided  by  BRL.  The  flat  base  case  was  selected  for  the  remaining 
calculations  to  reduce  geometry  related  complications  in  the  combustion  model 
development  This  grid  has  70  radial  points  from  the  centerline  to  the  projectile  base 
comer  and  120  axial  points  from  the  base  to  the  outflow  boundary.  A  modified  grid  for 
the  base  region  only  was  constructed  using  this  grid.  The  reacting  flow  calculations 
required  more  grid  points  in  both  the  radial  and  axial  directions  than  either  nonreacting 
or  air  injection  cases.  The  grid  utilized  for  the  base  calculations  has  100  points  radially 
from  the  centerline  to  the  base  comer,  and  69  points  fi’om  the  comer  to  the  outer 
boundary.  There  are  30  points  in  the  injection  region  which  extends  from  the  centerline 
to  a  radius  of  0.16  caliber.  In  the  axial  direction  there  are  46  points  upstream  of  the  base, 
and  150  points  from  the  base  to  the  outflow  boundary.  The  computational  grid  used  in 
the  base  flow  calculations  is  shown  in  Fig.  12.  The  near-wall  grid  resolution  in  the 
transverse  direction  (i.e.,  Ar/D)  was  2x  lO"* ,  which  provides  sublayer  resolutiotL  In  the 
axial  direction,  the  spacing  at  the  projectile  base  (i.e.,  AXg/D)  was  2x  lO"* ,  where  D  is  the 
maximum  projectile  diameter.  This  axial  spacing  is  about  7.5  times  smaller  than  that  in 
the  original  BRL  grid. 

Calculations  were  performed  for  the  M864  flat  base  projectile  for  the  following 
four  cases:  (a)  without  base  injection,  (b)  with  hot  air  injection,  (c)  with  H,  injection  and 
(d)  with  H  j  -CO  injection.  Case  (c)  used  the  global  reaction  model  of  Rogers  and 
Chinitz  (1983),  while  case  (d)  used  the  finite  rate  model  consisting  of  reactions  1  to  8  and 
10  to  12  in  Table  n.  For  case  (d)  the  injectant  species  mass  fraction  are  shown  in  Table 
I,  while  for  case  (c)  the  mass  fraction  of  H,  was  selected  to  give  the  same  fuel  heat  of 
combustion.  The  diluent  species  in  the  injectant  was  H^  O.  Also,  the  molecular  weight 
of  the  injected  H ,  was  artificially  increased  to  yield  the  correct  fuel  molecular  weight. 
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Figure  13  shows  the  favorable  comparison  between  the  SRA  and  BRL  base  pressure 
distributions  without  base  injection.  The  differences  between  the  grid  distributions  for 
the  two  calculations  were  noted  above,  and  is  expected  to  account  for  some  of  the 
observed  minor  differences  in  base  pressure. 

The  number  of  iterations  required  for  convergence  of  the  reacting  base  flow 
calculations  was  about  5500,  which  was  much  larger  than  that  needed  for  the  Jarrett 
validation  case.  This  is  due  to  the  complicated  recirculation  and  mixing  patterns  which 
occur  in  the  base  flow  region.  Also,  grid  resolution  in  the  shear  layers  and  reaction  zone 
will  influence  the  convergence  rate,  and  no  attempt  was  made  to  optimize  the  grid  for 
the  present  cases.  All  cases  with  injection  presented  here  have  an  injection  parameter, 
1=  0.0022,  where 


I 


mi 


(55) 


and  mj  is  the  injection  mass  flow  rate.  This  is  a  typical  value  for  I  at  projectile  launch 
Danberg  (1990).  The  stagnation  temperature  of  the  injected  gas  in  both  cases  (b)  and 
(d)  was  1533  K,  while  for  case  (c)  it  was  1755  K.  Although  the  higher  injection 
temperature  for  case  (c)  will  have  a  direct  impact  on  the  base  drag  coefficient,  the 
qualitative  features  of  the  flow  are  representative.  The  molecular  and  turbulent  Schmidt 
numbers  for  case  (c)  were  both  1.0,  while  for  case  (d)  they  were  0.7  and  1.0,  respectively. 
Reduction  of  the  molecular  Schmidt  number  for  case  (c)  was  not  carried  out;  however, 
based  on  calculations  for  H,  -CO  injection,  the  drag  coefficient  is  expected  to  decrease 
slightly  (about  7  percent).  In  reality,  the  molecular  diffusion  of  a  hydrogen-air  mixture  is 
not  correctly  represented  by  Pick’s  law  and  the  Schmidt  number  approximation.  In  a 
turbulent  flow  the  errors  encoimtered  tend  to  be  small  since  the  molecular  mixing  is 
often  much  less  than  the  turbulent  mixing.  In  the  present  case,  the  effect  of  the  Schmidt 
number  change  is  not  totally  negligible;  however,  a  greater  concern  is  the  turbulence 
prediction  in  the  base  region  since  the  mixing  of  reactants  clearly  will  have  an  important 
influence  on  the  combustion  process  and  the  resulting  base  drag. 

The  base  drag  coefficient,  Cdq,  is  defined  in  terms  of  the  maximum  projectile 
diameter  D  and  the  average  base  pressure  as 


j  [P*  -  PJdA 
^  base 


Ad  =  xD2/4 


(56) 
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The  present  calculations  yielded  values  of  Cqe  for  cases  (a),(b),  (c),  and  (d)  of  0.167, 
0.138, 0.118  and  0.089,  respectively.  Clearly,  there  is  significant  reduction  in  the  base 
drag  as  a  result  of  hot  air  injection  and  combustion.  The  results  shown  here  are 
consistent  with  the  predictions  of  Nietubicz  and  Sahu  (1988)  with  both  no  injection 
(Cdb=  0.167)  and  hot  air  injection  (Cdb=  0.136).  With  the  ninth  reaction  included  in 
the  calculation  the  drag  coefficient  obtained  with  -CO  combustion  was  also  about 
0.089,  which  indicates  the  relative  unimportance  of  that  reaction  in  the  present 
application.  The  base  pressure  distributions  for  these  four  cases  are  compared  in  Fig.  14. 
A  considerable  flattening  of  the  base  pressure  distribution  is  evident  as  a  result  of  both 
hot  gas  injection  and  combustion. 

Temperature  contours  and  velocity  vectors  for  the  four  cases  are  presented  in  Figs. 
15  (a,b,c,d)  and  16  (a,b,c,d),  respectively.  From  Figures  15  and  16  it  can  be  seen  that  the 
effect  of  hot  air  injection  is  confined  to  the  region  very  close  to  the  base.  The  H, 
combustion  case  (c)  shows  a  larger  effect,  while  the  -CO  combustion  case  (d)  shows  a 
much  more  extensive  region  of  influence.  In  addition,  the  stand-off  distance  of  the 
recirculation  region  increases  only  slightly  in  the  combustion  cases  (c)  and  (d).  The 
effect  on  the  base  pressure  may  be  explained  as  follows.  The  injection  of  gas  into  the 
base  region  will  move  the  recirculation  downstream  even  with  cold  injectant  This 
process  will  result  in  a  downstream  shift  in  the  wake  closure  location  (the  viscous  throat) 
and  a  somewhat  reduced  expansion  at  the  base  comer,  which  implies  an  increase  in  base 
pressure. 

With  hot  air  injection  the  temperature  increases  throughout  the  base  region  as  seen 
by  comparing  Figs.  16a  and  16b,  and  there  is  a  resulting  decrease  in  Mach  number. 

Hence  a  longer  distance  downstream  is  required  for  the  flow  to  reaccelerate  to  M  =  1 
and  the  viscous  throat  moves  further  downstream,  which  implies  a  weaker  comer 
expansion  and  higher  base  pressiue.  Strahle,  Hubbartt  and  Walterick  (1982)  explained 
this  behavior  in  terms  of  the  temperatme  increase  (Mach  number  decrease)  along  the 
stagnation  streamline. 

In  the  combustion  cases  the  wake  region  becomes  hotter  than  with  hot  air  injection. 
The  rear  stagnation  point  moves  even  further  downstream  and  the  comer  expansion 
becomes  even  weaker,  which  results  in  an  additional  increase  in  base  pressure.  With 
combustion  (case  c)  the  flame  zone  is  very  close  to  the  injection  plane,  and  the 
temp*  ratures  reached  are  higher  than  with  Hj-CXD  combustion  (case  d).  In  the  latter 
case  the  region  of  combustion  is  more  distributed,  and  temperatures  are  elevated  over 
much  more  of  the  projectile  wake  region.  The  effect  on  the  base  comer  expansion  and 
the  recompression  region  as  a  result  of  hot  air  injection  and  combustion  is  readily 
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observed  by  comparing  the  free  stream  temperature  contour  in  each  case.  These  are 
shown  in  Fig,  17  along  with  the  rear  stagnation  point  location  for  all  cases. 

53  Mesh  Refinement  Study  for  Reacting  Flow.  A  number  of  axial  grid 
refinement  calculations  were  performed  on  the  flat  base  case  with  combustion.  This  was 
done  to  assess  the  influence  of  axial  mesh  spacing  in  the  base  region  reaction  zone.  The 
importance  of  the  axial  spacing  at  the  projectile  base  and  near  the  downstream 
stagnation  point  was  also  identified.  The  original  case  considered  in  the  previous  section 
(case  d)  had  an  axial  spacing  axq/D  =  0.00Q2  at  the  base.  The  subsequent  grids  used  a 
relaxed  spacing  of  Axq/D  =  0.0015  at  the  base  to  permit  the  concentration  of  points  in 
the  reaction  zone  without  using  an  excessive  number  of  points.  The  latter  :>pacing  is  the 
same  as  that  in  the  original.BRL  grid,  which  has  been  used  by  BRL  personnel  in  base 
flow  calculations. 

Table  IV  indicates  the  axial  grid  spacing  in  three  relevant  locations  downstream  of 
the  base:  (1)  axq/D  is  the  spacing  at  the  projectile  base,  (2)  axi/D  is  the  spacing  in  the 
vicinity  of  the  first  stagnation  point  (closest  to  the  base),  and  (3)  axj/D  is  the  spacing  in 
the  vicinity  of  the  downstream  stagnation  point 

One  conclusion  which  can  be  reached  from  these  results  is  the  importance  of  the 
axial  spacing  at  the  base  upon  the  osdllatoiy  nature  of  the  base  pressure  convergence 
history.  The  drag  coefficients  reported  in  the  osdllatoiy  cases  (g,  h  and  i)  are  an  average 
taken  from  the  plot  of  versus  time  step  number.  With  this  in  mind  it  appears  that  a 
reduced  reaction  zone  spacing  characterized  by  Axj/D  is  not  critical  to  the  base  drag 
value;  however,  it  is  possible  that  overall  convergence  would  inq)rove  if  Ax^/D  was  also 
reduced  to  the  case  (d)  value.  Also,  the  cases  with  too  large  a  pacing  in  the  downstream 
stagnation  region  (ax^)  had  significant  convergence  problems. 

5.4  Two-Phase  Reacting  Flow  Case.  A  sample  two-phase  flow  calculation  has 
been  considered  with  injection  of  ammonium  perchlorate  (AP)  from  the  projectile  base. 
Actually,  the  partide  injection  was  initiated  slightly  downstream  of  the  gas  injection 
location  (at  xfD  =  5.925)  to  avoid  transporting  the  partides  through  the  Eulerian  fine 
grid  region  adjacent  to  the  base.  The  particle  mass  flux  was  assumed  to  be  one  per  cent 
of  the  injection  gas  mass  flux,  and  the  initial  partide  diameter  was  5.0  microns.  The 
effect  of  the  AP  burning  on  the  base  combustion  process  is  primarily  to  supply  additional 
O  3  in  the  base  region;  however,  the  relatively  low  AP  mass  flux  and  the  external 
freestream  air  supply  implies  that  the  effect  on  the  base  combustion  process  should  be 
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small.  This  was  observed  in  the  calculation  which  was  performed  by  restarting  the 
converged  gas  (Hj/CO)  combustion  case  and  running  an  additional  100  time  steps. 
Particles  of  diameter  less  than  1.0  micron  were  removed  from  the  calculation.  Figure  18 
shows  traces  of  representative  particles  injected  as  described  above.  These  traces  show  a 
fairly  slow  burnout  of  the  AP  particles  with  the  assumed  linear  regression  rate.  Note  that 
the  injection  gas  composition  was  not  modified  for  the  particle  injection  case.  If  a 
significant  fraction  of  the  propellant  AP  did  not  bum  in  the  combustion  chamber,  then 
the  equilibriim  composition  of  the  combustion  products  would  have  to  be  recalculated 
to  provide  new  gas  injection  boundary  conditions. 

6.  CXDNCLUDING  REMARKS 

Both  global  and  detailed  finite  rate  chemistry  models  have  been  incorporated  into 
an  existing  Navier-Stokes  computer  code  and  the  result  has  been  applied  to  the  problem 
of  combustion  in  the  base  region  behind  a  projectile  in  supersonic  flight  Prior  to 
considering  base  flow  and  combustion  cases,  the  code  was  calibrated  against  a  code  mn 
by  Ballistic  Research  Laboratory  and  when  similar  grid  and  turbulence  models  were 
used,  similar  results  were  obtained.  The  reacting  flow  capability  was  validated  through 
comparison  with  a  documented  supersonic  flow  coaxial  burner  experiment  When 
considering  the  effects  of  base  injection/combustion,  four  cases  were  simulated.  These 
were:  (a)  no  injection,  (b)  hot  air  injection,  (c)  H,  base  combustion  and  (d)  Hj/CO 
base  combustion.  The  predicted  base  pres.sure  increase  (base  drag  decrease)  is 
progressively  larger  with  hot  air  injection,  base  combustion,  and  -CO  combustion. 
The  detailed  computations  show  the  base  pressure  variation  to  be  intimately  related  to 
the  effect  on  the  base  comer  expansion  and  the  wake  recompression  region.  With  only 
injection  higher  flame  temperatures  are  reached  very  near  the  projectile  base  and 
lower  temperatures  occur  in  the  wake  region.  The  inclusion  of  CO  in  the  injected  fuel 
retards  the  reaction  process  and  lowers  the  peak  temperature  reached  near  the  projectile 
base.  However,  the  more  distributed  nature  of  the  base  combustion  yields  higher 
temperatures  throughout  the  wake  region  than  with  alone.  It  appears  reasonable  to 
include  both  and  CO  in  the  analysis  for  the  base  region  flow  physics  with  combustion. 

Preliminary  base  flow  calculations  employed  the  Jones  and  Launder  low  Reynolds 
number  k-c  turbulence  model,  while  subsequent  calculations  used  the  Baldwin-Lomax 
model  on  the  projectile  forebody  and  the  Chow  wake  mixing  model  downstream  of  the 
projectile  base.  The  k-e  model  calculations  showed  a  much  flatter  base  pressure 
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distribution  than  the  Chow  model,  although  this  may  be  a  result  of  performing  the  k-t 
simulation  on  a  grid  with  which  the  axial  resolution  next  to  the  base  adequate  for  a 
mixing  length  approach  but  too  coarse  for  a  viable  k-e  approach,  indicating  the  need  for 
considerably  higher  near  base  resolution  if  k-c  were  to  be  used.  Uncertainties  in  the 
base  region  turbulence  model  for  nonreacting  flows  appear  to  have  a  second  order  effect 
on  the  base  drag  coefficient  as  shown  by  both  present  computations  and  those  by  BRL 
personnel.  An  axial  mesh  refinement  study  for  the  reacting  base  region  flow 
demonstrated  the  importance  of  adequate  axial  resolution  in  the  vicinity  of  the 
downstream  stagnation  (wake  closure)  point  Additional  axial  grid  refinement  in  the 
combustion  region  over  that  in  the  baseline  case  did  not  significantly  change  the  results 
or  improve  the  convergence  rate  of  the  solution. 

A  validation  calculation  was  performed  for  a  coaxial  hydrogen-air  supersonic 
burner,  and  the  comparison  with  the  available  experimental  data  was  very  good.  This 
case  demonstrates  the  validity  of  the  reaction  set  which  is  a  subset  of  the  H,  -CO 
reaction  set  used  in  the  projectile  base  combustion  simulations.  Direct  validation  of  the 
base  combustion  models  developed  is  not  possible  without  additional  experimental  data 
for  base  combustion  flows.  Since  such  data  is  presently  unavailable,  the  validity  of  the 
computed  results  can  only  be  inferred  from  the  comparison  of  base  drag  coefficients  with 
those  implied  by  transonic  range  data  for  an  actual  projectile  firing. 

Further  work  should  be  performed  in  the  area  of  turbulence  -  combustion 
interaction  to  better  understand  the  impact  of  turbulence  unmixedness  on  the  H2  -CO 
combustion  process  in  the  base  region  flow.  Validation  of  any  models  would  require  a 
well  documented  base  flow  experiment  with  detailed  flow  field,  temperature  and  species 
data.  Although  a  significant  experimental  effort  would  be  required,  both  mean  and 
fluctuating  quantities  should  be  obtained  for  proper  model  validation. 
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Table  I.  M864  Propellant  Equilibrium  Species  Concentrations 
(Major  Species,  T  =  1533  K,  p  =  0.68  atm) 


Species  i 

CET86  results 

Equivalent  mixture 

mole  fraction 

Mwi 

mole  fraction 

mass  fraction 

CO 

0.249 

28.01 

0.265 

0.3402 

0.069 

44.01 

0.160 

0.3228 

HCl 

0.136 

36.46 

0.0 

0.0 

H, 

0.261 

2.016 

0.261 

0.0241 

H,0 

0.197 

18.015 

0.197 

0.1627 

N, 

0.085 

28.01 

0.117 

0.1502 

Sum 

0.997 

(Mw) 

21.81 

1.000 

1.000 

31 


Table  II.  Carbon  Monoxide  Oxidation  Mechanism  Including  HO, 
and  H,0,  Chemistry  (Westbrook,  et  al.  1977  and 
Gardiner  1984) 

Rate  constant:  kf  =  A  expC-Ejj/RuT)'*' 


Reaction 


A  b  Eg| 

(kJ/mole) 


1. 

H 

+ 

0* 

OH 

+ 

0 

1.2 

X 

10»T 

-0.91 

69.1 

2. 

H* 

0 

OH 

+ 

H 

1.5 

X 

10» 

2.0 

31.6 

3. 

0 

+  ' 

H.O 

OH 

+ 

OH 

1.5 

X 

io»» 

1.14 

72.2 

4. 

OH 

+ 

Ha 

«> 

H,0 

+ 

H 

1.0 

X 

10* 

1.6 

13.8 

5. 

0 

+ 

H 

+ 

M 

9 

OH 

+ 

M 

1.0 

X 

io»« 

0.0 

0.0 

6. 

O 

+ 

0 

+ 

M 

m 

0, 

+ 

M 

1.0 

X 

10*  » 

-1.0 

0.0 

7. 

H 

+ 

H 

+ 

M 

H, 

+ 

M 

9.7 

X 

10*  • 

-0.6 

0.0 

8. 

H,0 

+ 

M 

It 

H 

+ 

OH  + 

M 

1.6 

X 

10** 

0.0 

478.0 

9. 

Oj 

+ 

H, 

OH 

+ 

OH 

7.94 

X 

10*  < 

0.0 

187.0 

10. 

CO 

+ 

OH 

COa 

+ 

H 

4.4 

X 

10* 

1.5 

-3.1 

11. 

CO 

+ 

0 

+ 

M 

«t 

CO, 

+ 

M 

5.3 

X 

10*  » 

0.0 

-19.0 

12. 

CO 

+ 

Oa 

« 

CO, 

+ 

0 

2.5 

X 

10* » 

0.0 

200.0 

13. 

H 

+ 

0, 

+ 

M 

9 

HO, 

+ 

M 

2.0 

X 

10* » 

-0.8 

0.0 

14. 

0 

+ 

OH 

+ 

M 

9 

HO, 

+ 

M 

1.0 

X 

10** 

0.0 

0.0 

15. 

HO, 

+ 

0 

9 

0, 

+ 

OH 

2.0 

X 

10** 

0.0 

0.0 

16. 

H 

+ 

HO, 

9 

OH 

+ 

OH 

1.5 

X 

10** 

0.0 

4.2 

17. 

H 

+ 

HO, 

9 

H* 

+ 

Oa 

2.5 

X 

10** 

0.0 

2.9 

18. 

OH 

+ 

HO, 

9 

H,0 

+ 

Oa 

2.0 

X 

10** 

0.0 

0.0 

19. 

HO, 

+ 

HO, 

9 

HaO, 

+ 

Oa 

2.0 

X 

10** 

0.0 

0.0 

20. 

H,0, 

+ 

M 

9 

OH 

+ 

OH  + 

M 

1.2 

X 

10** 

0.0 

190.0 

21. 

H 

+ 

H,0 

2 

9 

HO, 

+ 

Ha 

1.7 

X 

10* » 

0.0 

15.7 

22. 

CO 

+ 

HO, 

9 

CO, 

+ 

OH 

1.5 

X 

10** 

0.0 

98.7 

23. 

H,0, 

+ 

OH 

9 

H,0 

+ 

HO, 

7.0 

X 

10** 

0.0 

6.0 

y  Set  B 


Set  A 


Dimensions  of  kf  are  [cm* /mole)  (®”1)  where  m  is  the 
reaction  order;  T  is  in  Kelvin.  Reverse  rate  constants  k], 
are  obtained  from  kf  and  the  equilibrium  constant  K^.. 
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Table  III.  Exit  Conditions  for  Jarrett  SSB  Coaxial  Streams 
(Jarrett,  et  al.  1988) 


Hydrogen 

Jet 

Vitiated 
Air  Jet 

Ambient  Air 

Mach  nuiaber 

1.0 

2.0 

0.0 

Temperature,  K 

545 

1180 

300 

Velocity,  m/s 

1772 

1390 

0.0 

Pressure, ,MPa 

Mass  Fractions 

0.112 

0.107 

0.101 

Y(H,) 

1.0 

0.0 

0.0 

0.0 

0.254 

0.232 

Y(N,) 

0.0 

0.572 

0.768 

Y(H,0) 

0.0 

0.174 

0.0 
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Table  IV.  Sununazy  Computed  Results  for  Projectile  Base 
Combustion 


Case 

Axial  Points 

AX3/D 

AX,/D 

AX^/D 

Cd3  and  Comments 

d 

196 

0.0002 

0.0136 

0.02 

0.089;  Nearly  monotonic 
convergence ; 

5500  steps. 

e 

msm 

0.0015 

0.007 

0.10 

-  Slow  convergence. 

f 

186 

0.0015 

0.003 

0.18 

-  Divergent. 

g 

186 

0.0015 

0.003 

0.19 

0.092  Erratic 

oscillatory 
convergence ; 

6000'‘'  steps. 

h 

186 

0.0015 

0.003 

0.05 

0.085  Oscillatory 
convergence ; 

OOOO"*"  steps. 

i 

196 

0.0015 

0.0035 

0.025 

0.092  Oscillatory 
convergence ; 

TOOO"*"  steps. 
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127  mm 


Janett  SSB  Experiment  (a)  Schematic  of  the  Apparatus 
(b)  Oomputational  Domain. 


Hgure  2.  101  X  101  Grid  for  Jairett  Supersonic  Coaxial  Burner  Simulation. 
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RADIAL  DISTANCE  (MM) 


Jarrett  SSB  Simulation  -  Temperature,  101  x  101  Grid. 


RADIAL  DISTANCE  (MM) 


JaiTctt  SSB  Simulation  -  Oo  Number  Density,  101  x  101 


RADIAL  DISTANCE  (MM) 


Jarrett  SSB  Simulation  -  N2  Number  Density,  101  x  101  Grid 


Figure  7.  Jarrett  SSB  Simulation  -  Temperature,  101  x  61  Grid. 


DIMENSIONS 


LENGTH  OF  PROJECTILE 

CALIBERS 

5.79 

NOSELB^TH 

CALIBERS 

3.42 

CYUNDER  LENGTH 

CALIBERS 

1.80 

BOATTAIL  LENGTH 

CALIBERS 

0.50 

BOATTAIL  ANGLE 

CEGPEES 

3.00 

Rgure  8.  Projectile  Schematic  (from  Danberg,  1990). 
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Surface  Pressure  (p/p^) 


0  1  2  3  4  5  6 

Axial  Distance  (calibers) 


Figure  10.  Forcbody  Surface  Pressure  Distribution  for  M864  Projectile 

with  Flat  Nose,  150  x  280  Grid  Symbols  from  BRL  Calculation- 
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Base  Pressure  (p/pj) 


l.d 

0.9 

0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.0  0.1  0.2  0.3  0.4  0.5 

Radial  Coordinate  (calibers) 

Figiue  11.  Base  Pressure  Distributions  for  Flat  Base  M864  Projectile. 

Present  Results  with  k-e  Turbulence  Model . 
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INFLOW  BOUNDARY 


CENTERLINE 


Figure  12.  169  x  196  Grid  for  M864  Flat  Base  Projectile  for  Region  Near  the  Base. 
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Base  Pressure  (p/p^) 


Base  Pressure  (p/p^) 


Radial  Coordinate  (calibers) 


Hguie  14.  Base  Pressure  Distributions  for  Cases  a-d  with 
Baldwin-Loiiiax/Qx)w  Turbulence  Model . 
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CENTERLINE 


Figure  15a.  Temperature  Contours  for  Case  (a):  Moo  =  2, 1  =  0.0,  Too  =  294  K,  T^  -  294 


Figure  15b.  Temperature  Contours  for  Case  (b):  Hot  Air  Injection,  Moo  =  2, 1  =  0.0022, 
To.-294K,Tw  =294K,Toinj“1533K. 
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Figure  15c.  Temperature  Contours  for  Case  (c):  H2  Injection,  Moo  =  2. 1  -  0.0022, 
Too  “  294  K,  T^  s=  294  K,  Tq  Jnj  =  1533  K. 


Hgure  15A  Temperature  Contours  for  Case  (d):  H2-CP  Injection,  Moo  =  2, 1  =  0.0022, 
Too  =  294  K,  T^  «  294  K,  Tq  inj  =  1533  K. 
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Figure  1 6b.  Velocity  Vectors  for  Case  (b):  Hot  Air  Injection,  Mo,  =  2, 1  =  0.0022, 
Too  =  294  K,  =  294  K,  To  inj  =  1533  K. 


Figure  16c.  Velocity  Vectors  for  Case  (c);  H2  Injection,  Mo©  =  2, 1  =  0.0022, 
To©  =  294  K,Tw  =  294  K,  Tq  i„j  =  1533  K. 


Figure  16d.  Velocity  Vectors  for  Case  (d):  H2-OO  Injection,  M©©  =  2, 1  =  0.0022, 
T©©  =  294K,Tw  =294K,Toi„j  =  1533K. 
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PROJECTILE 

BASE 


CASE  a 


CENTERUNE 


Figure  17.  Free  Stream  Temperature  Contoun  and  Rear  Stagnation  Points 
for  Cases  (a,  b,  c,  d). 
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Figure  1 8.  Representative  Particle  Traces  for  Two-Phase  Reacting  Flow. 
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UST  OF  SYMBOLS 


A  Area 

Ajy  Maximum  Projectile  Cross-Sectional  Area 

Af  Arrhenius  Constant  for  Reaction  r 

A  See  Eq.  (48) 


As 

Ajjase 

br 

C 


Cdb 

Ckieb 

Cpi(T) 

A 

<T)i 

Cs 

Qwk 

C, 

C, 

D 

D 

Djt 

Er 

F 

Fi 

Fi 

fi 

Gi 

H(«) 


Element  of  Time  Derivative  Matrix 

SeeEq.(49) 

Area  of  Projectile  Base 
Temperature  Constant  for  Reaction  r 
Defined  by  Eq.  (A-12) 

Factor  for  Baldwin-Lomax  Turbulence  Model 
Base  Drag  Coefficient 

Factor  for  Baldwin-Lomax  Turbulence  Model 

Frozen  Specific  Heat 

Specific  Heat  Per  Unit  Mass  of  Species  i 

Molar  Specific  Heat  of  Species  i 
Scaling  Parameter;  Eq.  (50) 

Factor  for  Baldwin-Lomax  Turbulence  Model 
Constant  in  Jones-Launder  Turbulence  Model 
Constant  in  Jones-Launder  Turbulence  Model 
Constant  in  Jones-Launder  Turbulence  Model 
Diffusion  Coefficient;  Maximum  Projectile  Diameter 
Nonlinear  Spatial  Difference  Operator  Matrix 

van  Driest  Damping  Factor 

Rate  of  Strain  Tensor 

Activation  Energy  of  Reaction  r 

Factor  for  Baldwin-Lomax  Turbulence  Model 

Conditioning  Factor;  Eq.  (SO) 

See  Eq.  (A-7) 

Mass  Fraction  of  Species  i  in  AP  Particle;  Eq.  (51) 
See  Eq.  (A-8) 

Time  Term  Vector 
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h 

ho 

hfi 

hi(T) 

I 

J 

Ji 

K 

k 

^>4- 

Kc^ 

t 

I 

mj 


Enthalpy  Per  Unit  Mass 
Stagnation  Enthalpy 

Heat  Formation  Per  Unit  Mass  of  Species  i 
Enthalpy  Per  Unit  Mass  of  Species  i 

Molar  Enthalpy  of  Species  i 
Injection  Parameter 
Jacobian 

Diffusive  Mass  Flux  at  Species  i 

Constant  in  Baldwin-Lomax  Turbulence  Model 

Turbulence  Kinetic  Energy 

Backward  Rate  Constant  for  Reaction  r 

Equilibrium  Constant  for  Reaction  r 

Forward  Rate  Constant  for  Reaction  r 

SeeEq.(48) 

Total  Number  of  Species;  Prandtl  Mixing  Length 

Injection  Mass  Flow  Rate;  Rate  of  Production  of  Species  i  Due  to 

Chemical  Reaction 


niy 

Ns 

n 

P 

Pi 

Po 

Pr 

PT"!* 

q 

qd 

Rcp 

Ri 

Rp 

R. 

r 

S 

S(«) 


Particle  Vaporization  Rate 

Total  Number  of  Species  in  the  System 

Power  for  Cartesian  and  Cylindrical  Polar  Coordinates 

Pressure 

See  Eq.  (A-9) 

Reference  Pressure  - 1  atm 
Molecular  (laminar)  Prandtl  Number 
Turbulent  Prandtl  Number 
Multicomponent  Energy  Flux 
Interdiffusional  Energy  Flux 
Particle  Reynold’s  Number;  Eq.  (54) 

Width  Parameter  for  Eggers  Turbulence  Model 
Particle  Radius 

Time  Rate  of  Change  of  Particle  Radius 

Universal  Gas  Constant 

Radius 

Defined  by  Eq.  (A-11) 

Source  Term  Vector 
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Sc 

Sop 

Sh 

T 

Tf 

To 

t 

U 

Ui 


Molar  Entropy  of  Species  i 
Molecular  (laminar)  Schmidt  Number 
Turbulent  Schmidt  Number 
Sherwood  Number;  Eq.  (53) 
Temperature 
Reference  Temperature 
Reference  Temperature;  298.15  K 
time 

Velocity  Vector 
Velocity  Component 


Ui'Uj' 

UdiEf 


Ut 

v’;,and*/"i, 

W 

Wi 

Xi 

[Xil 


Reynolds  Stress 

Difference  Between  Maximum  and  Minimum  Velocity  in  a  Shear 
Layer 

Wall  Shear  Velocity 

Stoichiometric  Coefficients  for  Reaction  r 
Defined  by  Eq.  (A-6) 

Molecular  Weight  of  Species  i 
Chemical  Symbol  of  Species  i 
Concentration  of  Species  i 


Xi 

Yi 

yj 

Yn 

y+ 

z 

’ll 

S 

€ 

fi 


K 

^eff 

f^cff 


Cylindrical  Polar  Coordinates 
Mass  Fraction  of  Species  i 
Computational  Coordinate  in  j-direction 
Distance  Normal  to  the  Wall 
Nondimensional  Normal  Distance 
Streamwise  Distance 
Defined  by  Eq.  (A-5) 

Defined  by  Eq.  (A-5) 

Local  Boundary  Layer  Thickness 
Dissipation  of  Turbulence  Kinetic  Energy 
Defined  by  Eq.  (A-5) 

Circumferential  Angle 
von  Karman  Constant 
Effective  Thermal  Conductivity 
Molecular  or  laminar  Viscosity 
Effective  Viscosity 
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^0 

ftj 


p 


T 


‘’’w 

♦ 


1«1 


i 


max 

min 


r 

s 

St 

w 

wake 


j 

n 


Molar  Chemical  Potential 
Turbulent  Viscosity 
Mass  Density 
Particle  Density 

Constant  in  Jones-Launder  Turbulence  Model 
Constant  in  Jones-Launder  Turbulence  Model 
Stress  Tensor 
Stress  Tensor 
Wall  Shear 

Equivalence  Ratio  for  the  Overall  Reaction  Process 
Dependent  Variable  Vector 

Magnitude  of  Vorticity 

Subcripts 

Species  Number;  Coordinate  Direction 

Maximum 

Minimum 

Reaction  Number 

Species  Number 

Stoichiometric 

Wall  Value 

Wake  Value 

Free  Stream 


Superscripts 

Computational  Direction 
Time  Step  Number 
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9.  APPENDIX  A 


Governing  Equations 

The  set  of  governing  partial  differential  equations  which  are  considered  are  the 
ensemble-averaged  Navier-Stokes  equations.  Before  these  equations  can  be  incorporated  into  a 
computer  code,  a  coordinate  system  must  be  chosen.  The  governing  equations  can  then  be  cast  in 
a  form  reflecting  the  choice  of  the  coordinate  ^tem.  The  coordinate  system  for  the  present 
calculations  must  be  capable  of  treating  general  configurations.  Therefore,  the  governing 
equations  written  in  a  cylindrical-polar  coordinate  system  are  transformed  with  a  general 
Jacobian  transformation  of  the  form 


yj  =  yj(Xi,  X,,  X3)  (A-l) 

where  (x ^ ,  x, ,  x,  )  —  (r,  6,  z)  are  the  original  (ylindrical  polar  coordinates.  The  velocity 
components  remain  the  components  (U  ^ ,  U, ,  U,  )  in  the  (x ^ ,  x, ,  x,  )  coordinate  directions, 
respectively.  The  new  independent  variables  yi  are  the  computational  coordinates  in  the 
transformed  system.  The  transformation  (A-l)  has  been  utilized  in  the  computer  code  both  with 
and  without  the  axisymmetric  flow  assumption. 


The  governing  equations  may  be  transformed  into  the  computational  space  using  the  chain 
rule  and,  after  some  manipulation,  written  in  the  following  v'ompact  form: 


+  JS  +  JC 

where  J  is  the  Jacobian  of  the  inverse  transformation  given  by 


(A-2) 


a(ys  y’,  y’) 


(A-3) 


and 
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(A-4) 


yj  . 

^  ,t  at 


Further,  the  coefficients  7i,  f  i  are  given  by 


ITS  =  1 


(A-5) 


and  m  s  1  for  all  equations  except  the  X2-<lirection  momentum  equation,  for  which  m  »  2.  The 
vector  variables  used  in  Eq.  (A-2)  are  defined  as 


(A-6) 


pU,Ui 

pU,Ui 

pU,Ui 

pUi 

PhoUi 

pkUi 

pcUi 

pY^Ui 


(A-7) 
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^  ^11 
r^nr. 


-r^qi 


"(55  ■"  ^]’'A.i 


where  n  =  1  for  i  =  1  and  n  =  0  for  i  =  2, 3. 


(A-8) 


Pi  = 


(A-9) 


Note  that  the  velocity  components  (U ^ ,  U, ,  U,  )  are  the  (ylindrical-polar  velocity  components, 
and  T  nj  is  the  stress  tensor  written  in  cylindrical-polar  coordinates.  The  molecular  and  turbulent 
stress  tensor  may  be  written  as 


Tij  =  2/ieff  eij  -  -fieff 


(A-10) 


where  e^j  are  components  of  the  rate  of  strain  tensor  in  (^lindrical-polar  coordinates.  The 
multicomponent  energy  flux,  qj,  is  given  by  Eq.  (7)  of  the  main  text. 

Finally,  the  vector  S  contains  source  terms  and  certain  differential  terms  which  do  not 
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conform  to  the  basic  structure  of  Eq.  (A-2),  and  the  vector  C  contains  the  additional  curvature 
terms  due  to  the  cylindrical-polar  coordinate  ^tem. 


nV 

0 

0 


8  = 


^  +  v.(r.O) 


2^T  C  ®ij  ®ij 


^  —  2pv 


2/i/iT 


Cl  ^  /^T  [  2eij  eij  ]  +  [v^d]  - 


CjP  ^ 


L  ®k 


(A-11) 


I  pu/  -  I  T,, 


-  r 


0 


c 


0 


0 

0 

0 

0 


The  chemical  source  term,  mjp  is  determined  from  Eqs.  (17)  and  (22), 


(A-12) 


d[XK) 
®k  “  '^k  dF" 


(A-13) 
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